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ABSTRACT 

We have developed an improved recombination calculation of H, He I, and He II in the early Universe 
which involves a line-by-line treatment of each atomic level. We find two major differences compared 
with previous calculations. Firstly, the ionization fraction Xc is approximately 10% smaller for redshifts 
<^ 800, due to non-equilibrium processes in the excited states of H. Secondly, He I recombination is 
much slower than previously thought, and is delayed until just before H recombines. We describe the 
basic physics behind the new results and present a simple way to reproduce our calculation. This should 
enable fast computation of the ionization history (and quantities such as the power spectrum of CMB 
anisotropics which depend on it) for arbitrary cosmologies, without the need to consider the hundreds 
of atomic levels used in our complete model. 

Subject headings: cosmology: theory — atomic processes — early universe — cosmic microwave 
background 



1. INTRODUCTION 

We have recently revisited the calculation of the re- 
combination epoch of the early Universe by making as 
few approximations as possible and trying to retain full 
accuracy (Seager, Sasselov, & Scott 1999, hereafter Pa- 
per I). We were motivated by the potential to measure 
(with the MAP and Planck satellites) CMB anisotropics 
at the roughly 1% level over a wide range of angular scales, 
and indications (Hu et al. 1995, hereafter HSSW) that the 
existing solution for hydrogen and helium recombination 
has uncertainties at that level. Indeed, we uncovered a 
number of minor improvements, as well as two more major 
effects which change the ionization history in a significant 
way. 

Detailed understanding of the recombination process is 
crucial for modeling the power spectrum of CMB aniso- 
tropics. Since the seminal work of the late 1960s (Peebles 
1968, Zel'dovich et al. 1968), several refinements have been 
introduced (see discussion in HSSW and Paper I), but lit- 
tle has changed. Modern codes for evolving the ionization 
fraction = ric/nn (where rie is the number density of 
electrons and tt-h is the total number density of H nuclei) 
have been based almost entirely on the single differential 
equation introduced 30 years ago, with a more accurate re- 
combination coefficient, but no other basic improvement. 

With today's computing power, there is not the neces- 
sity for making sweeping approximations that existed 30 
years ago. We believe our work represents the most ac- 
curate picture to date of how exactly the Universe as a 
whole became neutral. In this Letter we summarize what 
is new in the physics involved and present an approximate 
treatment of our set of equations which accurately repro- 
duces our complete calculation. We discuss in detail the 
physical basis for each approximation and pay attention 
to the limits of validity. The paper is supplemented by 
a computer code (reef ast) which can be used to do the 
calculation, and in conjunction with, e.g. cmbf ast by Zal- 



darriaga, Spergel & Seljak (1997), to compute accurate 
CMB power spectra for different cosmologies. 

2. HOW THE UNIVERSE BECAME NEUTRAL 
2.1. Our Multi-level Calculation 

In the canonical Hot Big Bang picture, the recombi- 
nation epoch is when the Universe became cool enough 
for protons to capture electrons and form neutral hydro- 
gen. This recombination process was not instantaneous, 
because the electrons, captured into different atomic en- 
ergy levels, could not cascade instantaneously down to the 
ground state. The electrons were impeded because of fast 
reionizations out of excited states due to the huge reser- 
voir of low energy photons, and because of the high op- 
tical depth of the Lyman lines and continuum transitions 
to ground state. Any Lyman line or continuum transi- 
tion to the ground state emitted a photon with energy 
where there were few blackbody photons, which immedi- 
ately photoexcited or photoionized a neighbouring atom in 
the ground state. Atoms reached the ground state through 
either cosmological redshifting of the Ly a line photons, or 
by the 2s-ls two photon process. Because these rates 
from n=2 to the ground state were much slower than the 
net recombination rate to n=2, a "bottleneck" occured 
which slowed down the entire recombination process. The 
Universe expanded and cooled faster than recombination 
could be completed, and a small fraction of free electrons 
and protons remained. This fraction, during and after 
recombination, affects the CMB anisotropics through the 
precise shape of the thickness of the photon last scattering 
surface (i.e. the visibility function). 

The 'standard' methodology considers an 'effective 3- 
level atom' with a ground state, first excited state (n=2) 
and continuum, with the states n > 2 represented by a 
recombination coefficient. A single ordinary differential 
equation can then be derived to describe the ionization 
fraction, equation (Q) (See Peebles, 1968, 1993). Many 
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assumptions go into this derivation, including: that H ex- 
cited states are in equilibrium with the radiation; that 
stimulated deexcitation is negligible for the Lya transi- 
tion; that a simple recombination coefhcient can be used; 
that every net recombination results in a ground state 
atom, so that the ground state number density ni — 
"■H — n-p (where tih is the total hydrogen number density, 
ionized or neutral, while Up is the proton number density, 
i.e. ionized H); that the Lya redshifting can be dealt with 
using a simple escape probability; that collisional processes 
are negligible; and that He can be ignored. 

Our new methodology, made possible with modern com- 
puting power, is to calculate recombination with as few 
approximations as possible. Instead of evolving a sin- 
gle ODE for dxc/dt for H recombination, we evolve one 
ODE for each of 300 atomic energy level populations, as 
well as one for electrons and one for the matter tempera- 
ture. With this level- by- level treatment, we include thou- 
sands of bound-bound transitions which correctly couple 
the ODEs, and we calculate recombination to each atomic 
level as it evolves with redshift which eliminates the need 
for an effective recombination coefficient [a in §3). We 
also included other effects (which turned out to be neg- 
ligible), such as collisional transitions, feedback of distor- 
tions to the radiation field, complete heating and cooling 
terms, and H chemistry. See Paper I for full details on this 
method. Our new method tests the approximations used 
in deriving the standard ODE, and we find all of them 
to be valid, except for the equilibrium assumption. The 



causes and consequences of this are discussed in §2.2.1 



and the effect can be approximated using the single ODE 
approach, with a small modification that artificially speeds 
up the recombination at low z, as discussed in 

With larger ionization potentials. He I and He II re- 
combined before H. He II recombination shows no de- 
viation from previous calculations, in that we find only 
small differences compared with Saha equilibrium. How- 
ever, He I recombination is important for the CMB an- 
isotropics (HSSW), and to some extent for the chemistry 
in the early Universe (Standi et al. 1999), and requires 
closer scrutiny. Our code evolves He I and He II in the 
same level- by- level method describe above, and simulta- 
neously with H. Previous calculations for He I used the 
Saha equilibrium method (e.g. HSSW) or an effective 3- 
level method (Matsuda, Sato & Takeda 1971). We find 
that the excited states of He I remain in equilibrium with 
the radiation (unlike H), but that accurate treatment of 
He I results in a much delayed recombination. This can 
also be reproduced using a simple ODE for He I. 

2.2. What is New in the Physics 

2.2.1. Hydrogen 

Neutral H production is slightly faster when using our 
multi-level calculation, compared to the standard calcula- 
tion (updated with the m ost recent parameters (HSSW)), 



as shown in Figure 2.2.1 



For the standard equilibrium case, the net bound-bound 
rates are zero, and this is an implicit assumption in deriv- 
ing equation (|^). We find that at z ^ 1000, the net bound- 
bound rates become different from zero, because at low 
temperatures, the cool blackbody radiation field means 
there are few photons for photoexcitation of high energy 



transitions (e.g. 70-10, 50-4 etc.). In this case sponta- 
neous deexcitation dominates, causing a faster downward 
cascade to the n=2 state. In other words, once an elec- 
tron is captured to, say, n=70, it can cascade down to 
the n=2 state faster than in the equilibrium case, because 
few photons are around to photoexcite it. In addition, the 
faster downward cascade rate is faster than the photoion- 
ization rate from the upper state, and one might view this 
as radiative decay stealing some of the depopulation "flux" 
from photoionization. Both the faster downward cascade 
and the lower photoionization rate contribute to the faster 
net recombination rate. 



500 



1000 
redshift 



effective 5- 
300-level 



1500 



2000 



Fig. 1: Multi-level hydrogen recombination for the standard 
CDM parameters fitot = 1.0, = 0.05, Hq = 50, Fp = 
0.24, To = 2.728 K. The two lines are the complete calculation 
of Paper I (dotted) and the standard (effective 3-level atom) 
calculation (dashed). 

An imbalance develops as the background radiation con- 
tinues to cool while the entire Lyman series remains op- 
tically thick (the "bottleneck" gets worse). A Boltzmann 
distribution relative to n=2 is no longer sustainable — the 
excited states become progressively overpopulated. (Note 
that this is not a population inversion.) The radiation 
field is cool, but strong, and able to keep neighbouring 
states well coupled, as well as to hold the highest Ry- 
dberg states in equilibrium with each other. However, in 
comparison with the standard equilibrium capture-cascade 
calculation for a, this unusual situation leads us to higher 
effective recombination rates to the majority of excited 
states, without increasing photoionization proportionally. 
This results in a higher net rate of production of neutral 
hydrogen atoms. 

Because the net effect of our new H calculation is a faster 
recombination (a lower freeze out ionization fraction), our 
results can be reproduced by artificially speeding up re- 
combination in the standard calculation, simply by mul- 
tiplying the recombination and ionization coefficient by a 
'fudge factor' F. 

2.2.2. Helium 

Our new multi- level calculation shows a significantly de- 
layed He I recombination compared to previous calcula- 
tions. In fact, for our low Qb models He I recombination 
is still finishing at the very beginning of H recombination. 
Figure 2.2.2 shows the ionization fraction Xe through He II, 
He I and H recombination, plotted against the standard H 
calculation and He II and He I recombination calculated 
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using the Saha ionization equilibrium equation. 
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Fig. 2: Helium and hydrogen recombination for the standard 

CDM parameters with Yp = 0.24 and To = 2.728 K. The first 
step from right to left is recombination of He HI to Ho H, the 
second step is He H to He I. The dashed line is our new He 
I calculation, while the dotted line is the He Saha equilibrium 
recombination and the H standard calculation. 

The recombination of He I is slow for the same rea- 
sons that H recombination is, namely because of the op- 
tically thick n^p-l^s transitions which make cascades to 
the ground state slow, and the exclusion of recombinations 
to the ground state. In other words He I follows a Case 
B recombination with an inhibition factor. Because the 
"bottleneck" at n=2 largely controls recombination, it is 
not surprising that He I and H recombination occur at a 
similar redshift; the ionization energy of n=2 is similar in 
both. 

The physics of He I recombination is different than that 
for H recombination because of its different atomic struc- 
ture. More specifically, the high excited states of He I arc 
much more strongly packed towards the continuum com- 
pared to those of H: the energy difference between the 3p 
levels and the continuum is 1.6 eV for He I versus 1.5 eV 
for H, compared to 24.6 cV versus 13.6 cV for the groimd 
state. Because of this, the radiation field has a stronger ef- 
fect on the excited states of He I than those of H. This has 
3 main consequences: 1) The strong radiation field during 
He I recombination keeps the ratio of photoionization rate 
to downward cascade rate higher than the H case, resulting 
in a slower recombination than H. 2) The strong radiation 
field also causes the triplet states (e.g. n^p-n^s) to be vir- 
tually unpopulated. The triplets have only forbidden (i.e. 
very slow) transitions with the singlets, and a metastable 
"ground state" 4.77 eV from the continuum. 3) The radia- 
tion field is strong enough that the excited states remain in 
equilibrium with the radiation throughout recombination. 

There is one possibility to speed up the slow He I re- 
combination. It is the existence of some neutral H, which 
could 'steal' He I resonance line photons, invalidate the ef- 
fective Case B by removing the "bottleneck" at n=2, and 
make it a Saha equilibrium recombination instead. How- 
ever, our detailed calculation shows that the rate for this 
process never comes close to being significant. 

Because He I recombination simply follows a Case B 
recombination and unlike for H the excited singlet states 
remain in equilibrium throughout recombination, it can be 
reproduced using the effective 3-level single ODE (shown 



below). Because the triplet excited states remain unpop- 
ulated during recombination, the effective 3-level atom is 
the ground state l^s, the first singlet excited state, and 
the continuum. 

3. THE APPROXIMATE SET OF EQUATIONS 

Here we present a set of equations which allow a simple 

way to reproduce our new recombination results from the 
multi-level code. Though approximate, the set provided 
below includes more of the recombination physics than 
the 'standard' calculation, and reproduces approximately 
(via parametrizations) the departures from equilibrium in 
H, and the slow He I recombination found in Paper I. 

However, a word of caution is necessary - the set of 
equations below have a range of validity which corresponds 
roughly to the choice of typical cosmologies discussed in 
Paper I. Attempts to calculate recombination in more ex- 
treme cosmologies, or with the addition of extra physics, 
will probably require evolution of the full set of equations 
of the multi-level code. In addition, our approximations 
here are designed for optimal use with the CMB aniso- 
tropies, and do not fit as well the range below z < 300, 
where molecular formation becomes important. Therefore 
detailed study of the chemistry in the early universe will 
probably also require evolution of the full set of equations 
of the multi-level code. 

The approximate set of equations below — two ODEs 
for the ionization fractions of H and He I, and an equation 
for the matter temperature — should be solved simulta- 
neously. We recommend leaving out He II recombination 
entirely, since it has no effect on the power spectriim of 
CMB anisotropics. Otherwise Saha equilibrium is an ade- 
quate approximation. The two ODEs derive from consid- 
eration of detailed balance in the effective 3-level atoms of 
H and He I (see Peebles, 1968, 1993). The total recom- 
bination coefficients in each case have been parametrized 
to reproduce our multi-level results. The equation of the 
total rate of change of the matter temperature is a trun- 
cated version of Paper I eq. (69), including adiabatic and 
Compton cooling terms. The set is as follows: 
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dz 3i7(z)(l + z)r7ioC 1 + /hc + a^c (1 + -2) 

. ^ . . (5) 

Equation (0) is for the matter temperature, which we 

recommend to use in the entire calculation above, due to 
the small but important effect resulting from the differ- 
ence between Tm and Tr at low redshift. Below we list 
all constants and parameters which appear in the above 
equations. 

The constants are Boltzmann's constant fc, Planck's con- 
stant /i, the speed of light c, Thomson scattering cross 
section ctt, the electron mass mo, and the radiation con- 
stant flR. The three independent variables are proton frac- 
tion x-p = fip/nu, singly ionized helium fraction xhoII — 
n^ieii/n-a, and matter temperature Tm- The dependent 
variable is the electron fraction Xc — nc/nn — Xp + xhcII- 
Here n refers to number density, and tih is the total hy- 
drogen number density. 

Turning to the atomic data, the H Lya rest wave- 
length is Ah2p = 121.5682 nm. The H 2s-ls frequency 
t^H2s = c/Ah2p, is close enough to Lya that the same av- 
eraged wavelength value can be used. The He I 2^p-l^s 
wavelength is X^cu^p = 58.4334 nm. Note that the He I 
2^s-l^s frequency i^hc/2s — c/60.1404nm. Unlike for H, 
the separation of He I 2^p and 2^s is large enough that 
AHe2ip and Ahc21s must be distinguished, hence the extra 
exponential term in equation (|^) over equation (|^), with 
t'He/2ip2is = ^'He/2ip " t'He/2is = t'ps (whcrc wc define this 
last notation to save space). The H 2s-ls two photon 
rate is Ah = 8.22458 s'^ (Goldman, 1989), while the He I 
2^s-l^s two photon rate is Ahc = 51.3 s~^ (Drake et al. 
1969). 

Further, an is the Case B recombination coefficient 
for H, from Hummer (1994), and fitted by Pequignot et 
al. (1991), with a = 4.309, b = -0.6166, c = 0.6703, 
d = 0.5300 and t = Tm/10'* K. The factor F is 1.14, and al- 
lows equation (|^) to agree with our multi-level calculation 
by speeding up recombination. Note that it also enters into 
the coefficient /? via a as described below, ane is the Case 
B He I recombination coefficient for singlets, from Hum- 
mer and Storey (1998). The parameters are q = 10"^^'^"''*, 
p = 0.711, Ti = lO^-^-^K, and T2 fixed arbitrarily at 
3K. This fit is good to < 0.1% over the relevant tempera- 
ture range (4,000-10,000 K), and still fairly accurate over a 
much wider range of temperatures. The P are photoioniza- 
tion coefficients and are calculated from the recombination 
coefficients by: /3 — a{2nmck Tm/^^)"^^^ exp(— /ij/2s/fcTM)- 
Here 1/23 and a are different for H and He I. Note that Tm 
and J^2s are used here, and that incorrectly using Tr or i'2p 
will cause a small but important difference for high baryon 
models. 

The cosmological parameters are the redshift z, Hubble 
factor H(z), the cosmological redshifting of H Lya photons 
Kn = /(87r_ff(z)), the cosmological redshifting of He I 

2^p-l^s photons i^Hci = A^jjj /{8ttH{z)), and the radi- 

ation temperature Tr = To(l + z). The primordial He 
abundance was taken to be Ip = 0.24 by mass (Schramm 
& Turner, 1998), and the present-day CMB temperature 
To to be 2.728 K (the central value determined by the FI- 
RAS experiment, Fixsen et al. 1996). 



A word of caution about the numerical computation. 
The equations are generally stiff, and so there are two 
approaches which can be taken. One is to use an integra- 
tor appropriate for stiff sets of equations, the other is to 
use Saha equilibrium values of to carry the integrations 
through the stiffest regimes (the beginning of each recom- 
bination epoch). The former case is not much slower, and 
certainly will work for arbitrary cosmology. For the latter 
case some experimentation may be necessary in order to 
efficiently choose the redshift to start and finish the Saha 
approximation for particular cosmologies, and then it may 
be faster. 

For He II it is sufficient to use the relevant Saha equa- 
tion: 

(xe - 1 - /He)a:c ^ {2nm^kT)^/^ 

l + 2fiie-Xe h^TlH 

For the non-stiff integrator approach, this can be switched 
on at say z > 5,000. Then there will be a period when 
He is all singly ionized, and after say z = 3500 the Saha 
equation for He I needs to be switched on: 

(^c - l)^e ^ ^ i2nm,kTf/' ^_^^^^,,^_ 
1 + /hc - h^nn 

This can be used until say 1% of the He I has recombined, 
at which point the coupled ODEs can be solved using a 
routine such as dverk (Hull et al., 1976). It is also worth 
fixing the hydrogen to be fully ionized until some redshift, 
then using Saha for H to get recombination started. It 
may also save some integration time if He recombination 
is switched off once the fraction of neutral He has fallen 
below some small number. 

We found that in fact it was not necessary to evolve the 
H and He equations simultaneously, at least for cosmolo- 
gies which we checked in detail. It appears to be sufficient 
to evolve each separately, even in cases where there is a 
small overlap in the recombination epochs. However, there 
is little computational expense in solving the 3 ODEs si- 
multaneously, and so that is what we recommend. 

4. CONCLUSION 

We have presented the basic physics behind our im- 
proved recombination calculation which shows a signifi- 
cantly delayed He I recombination and a 10% lower resid- 
ual Xe at freezeout compared to previous calculations. We 
present a set of equations to reproduce our work, which 
are modified versions of those previously used because our 
new, detailed calculation agrees very well with the results 
of the standard calculation. This underscores the tremen- 
dous achievement of Peebles, Zel'dovich and colleagues in 
so fully understanding cosmic recombination 30 years ago. 
However, the goal of modern cosmology is to determine 
the cosmological parameters to an unprecedented level of 
precision, and in order to do so it is now necessary to un- 
derstand very basic things, like recombination, much more 
accurately. 
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